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Abstract. 

Recent angle-resolved photoemission spectroscopy (ARPES) has identified that a 
finite-range Frohlich electron-phonon interaction (EPI) with c-axis polarized optical 
phonons is important in cuprate superconductors, in agreement with an earlier proposal 
by Alexandrov and Kornilovitch. The estimated unscreened EPI is so strong that 
it could easily transform doped holes into mobile lattice bipolarons in narrow-band 
Mott insulators such as cuprates. Applying a continuous-time quantum Monte-Carlo 
algorithm (CTQMC) we compute the total energy, effective mass, pair radius, number 
of phonons and isotope exponent of lattice bipolarons in the region of parameters 
where any approximation might fail taking into account the Coulomb repulsion and 
the finite-range EPL The effects of modifying the interaction range and different lattice 
geometries are discussed with regards to analytical strong-coupling/non-adiabatic 
results. We demonstrate that bipolarons can be simultaneously small and light, 
provided suitable conditions on the electron-phonon and electron-electron interaction 
are satisfied. Such light small bipolarons are a necessary precursor to high-temperature 
Bose-Einstein condensation in solids. The light bipolaron mass is shown to be universal 
in systems made of triangular plaquettes, due to a novel crab-like motion. Another 
surprising result is that the triplet-singlet exchange energy is of the first order in 
the hopping integral and triplet bipolarons are heavier than singlets in certain lattice 
structures at variance with intuitive expectations. Finally, we identify a range of 
lattices where superlight small bipolarons may be formed, and give estimates for their 
masses in the anti-adiabatic approximation. 
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1. Introduction 

A growing number of observations point to the possibility that high-Tc cuprate super- 
conductors may not be conventional Bardeen-Cooper-Schrieffer (BCS) superconductors 
[l], but rather derive from the Bose- Einstein condensation (BEC) of real-space pairs, as 
proposed by Mott and others [21 |3l HI [5] . A possible fundamental origin of such strong 
departure of the cuprates from conventional BCS behaviour is the unscreened (Frohlich) 
EPI with a polaron shift, Ep of the order of 1 eV (La2Cu04, Ep ^ 0.65eV) [6], routinely 
neglected in the Hubbard U and t — J models |[7j. This interaction with c— axis polar- 
ized optical phonons is virtually unscreened because the upper limit for the out-of-plane 
plasmon frequency (< 200 cm~^ [8]) is well below the characteristic frequency of opti- 
cal phonons, uj ~ 400 - 1000 cm Since screening is poor, the magnetic interaction 
remains small compared with the Frohlich EPI at any doping of cuprates. In order to 
generate a convincing theory of high-temperature superconductivity, one must treat the 
Coulomb repulsion and unscreened EPI on an equal footing. When both interactions are 
strong compared with the kinetic energy of carriers, the so-called "Coulomb-Frohlich" 
model (CFM) predicts a ground state in the form of mobile, preformed, inter-site pairs 
dressed by lattice deformations (i.e intersite bipolarons) [HllHlll]. 

The most compelling evidence for (bi)polaronic carriers in novel superconductors 
is the discovery of a substantial isotope effect on the carrier mass jlOj predicted by the 
bipolaron theory [11]. Recent high resolution ARPES [121 [13] provides another piece 
of evidence for a strong EPI in cuprates between electrons and c-axis-polarised optical 
phonons [13]. These, as well as recent tunnelling [14J, earlier optical [15] and neutron 
scattering [16] experiments unambiguously show that lattice vibrations play a significant 
though unconventional role in high temperature superconductors. 

Remarkably, earlier path-integral studies of large bipolarons in the continuous 
limit [IT]) led to a double surprise: (a) The large bipolaron is only stable in a very 
limited sector of the parameter space (Coulomb repulsion versus the Frohlich coupling 
constant) (b) Most traditional "Frohlich polaron" materials (alkali halides and the like) 
lie completely outside (and "far" from) this bipolaron stability sector, but several high- 
Tc superconductors lie very close and even inside this rather restricted area of stability 
in the parameter space. 

When the strong Frohlich EPI operates together with a shorter range deformation 
potential and molecular-type (e.g. Jahn- Teller) EPIs, it readily overcomes the Coulomb 
repulsion at short distances of about the lattice constant, so that large (continuous) 
bipolarons become local (lattice) bipolarons in narrow bands [4]. Even at significant 
doping local pairs are not overlapped, so that a high critical temperature for Bose- 
Einstein condensation (BEC) could be achieved, if they are sufficiently mobile. Analysis 
of the site-local Holstein-Hubbard model has indicated that in order for the Coulomb 
repulsion (Hubbard U) to be overcome by the induced attractive force between electrons, 
EPI must be so large that the polaron (and bipolaron) masses must be huge, rendering 
the transition temperature minuscule |18j. All is not lost, however, since the Holstein 
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interaction is the extreme short-range hmit of a finite range EPI. Using the finite-range 
EPI it is possible for electrons to pair between sites [6l |9] without requiring the electron- 
phonon induced attraction to be larger than the Hubbard U. Moreover, the individual 
polarons are significantly lighter, so the mass of the pair has potential to be orders of 
magnitude smaller than in the Holstein case. 

To put these arguments on a solid microscopic ground we simulate the CFM 
Hamiltonian on a lattice using an advanced QMC technique for bipolarons and compare 
with analytic results in the strong coupling and anti-adiabatic limits. First, we introduce 
the model. 



2. Coulomb- Prohlich model 



The Hamiltonian for the CFM is written as 

(nn')(T nn'aa' 
m m nmcr 

Each ion has a displacement ^m- Sites labels are n and m for electrons and ions 
respectively, c annihilate electrons. The phonons are Einstein oscillators with frequency 
uj and mass M. (nn') denote pairs of nearest neighbours, and = —ihd/d^-cn ion 
momentum operators. The instantaneous interaction l^(n, n') has on-site repulsion U 
and nearest neighbour interaction V (if the electron-phonon coupling term is set to zero, 
one obtains the simple UV model). The force function is of the screened Frohlich type, 

. , , K / |m — n| \ , . 

= [(m-n)^+i]3/^ ) 

(k is a constant) [20]. We will also use a slightly different notation for the electron- 
phonon interaction term here, H^i-ph = Xlijo- fl'«i4crCio-('^j + dj) where gij is a 
dimensionless interaction proportional to the force, and d}^ create phonons at site j. 
We set h=l. 

Such a model has a remarkable property. Unlike the site local Holstein model, 
there is attraction (and potentially pairing) even in the presence of very strong on- 
site Coulomb repulsion. The model is justified in the presence of alternating planes of 
itinerant electrons and ions, where there is strong screening along the c-axis. 

There have been a number of studies discussing the masses of polarons and 
bipolarons with long range interaction |T9l [2T], [22]. The polaron formed from the 
long-range Frohlich interaction proposed in [6] has been simulated in reference [9], 
demonstrating that the polaron mass may be significantly lighter than its Holstein 
counterpart. This is due to the nature of polarons in the Holstein case, which 
may be demonstrated nicely by examining the Lang-Firsov transformation. In that 
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Figure 1. (a) Staggered vs (b) rectangular ladders. Ions are placed one lattice spacing 
above a ladder of electrons, with one ion per site. The ions are permitted to vibrate 
in the z-direction only. Electrons inhabit one leg each, with no hopping between 
the legs. In the strong coupling limit, there are significant geometrical differences. 
On the staggered ladder, two degenerate near neighbour pairs (A and A') can form, 
which allows the polaron to scuttle in a crab-like manner with mass proportional 
to the polaron hopping. Alternatively, on the rectangular ladder there is one near- 
neighbour state, and in order to move, the pair must break to state A in order to change 
configuration from state B to B'. Such a state propagates by waddling awkwardly, and 
has mass proportional to polaron hopping squared. 



transformation, the operators in the Hamiltonian are replaced in the following way, 



thus hopping processes in the Holstein polaron (where gij = g6ij) take place by a 
complete relaxation of the lattice on the initial site, a hop, and then a distortion on 
the target site. With a longer range interaction, the lattice is pre-distorted before the 
particle moves, leading to a much smoother process with a lower intermediate energy 
state. We have recently determined that long-range interactions lead to a reduction 
of the importance of geometry on the properties of the bipolaron, especially the mass, 
leading to very similar results on triangular and square lattice [21]. We will discuss the 
crossover between Frohlich and Holstein polarons later in this article. 

The Lang-Firsov transformation is an exact canonical transformation, and leads 
to a transformed Hamiltonian with a new transformed wavefunction \'^)lf = e~'^|\E'). 
It is most instructive to consider the transformation of the atomic Hamiltonian and 
the transformation of the hopping terms separately, since typically the Lang-Firsov 
transformation is the starting point for a series of perturbative analyses. 
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2.1. Transforming the atomic Hamiltonian 

When the hopping term is set to zero, the phonon portion of the CFM is written as 
follows: 



Hat = ^ 



(N.B. The index i is now taken to contain a spin and a site index). Applying the 
Lang-Firsov canonical transformation 0, 



Hat= - X] dij^iid] + + 2 ^ gi'jrii'] 



(6) 



ll' j j 

where we reordered the summation, and noted that gij = fij/uV2Muj. This shows the 
remarkable property of the Lang-Firsov transformation, since the electron and phonon 
subsystems in the atomic limit are now completely decoupled. 

At this stage, it is convenient to introduce the following function, 

<|.Ar[r(r), r(r')] = /m[r(r)]/„.+Ar[r(r')] (8) 

m 

where the reason for adding an additional translation, Ar, to the phonon sub-system will 
become apparent when the action is introduced in the next section. For the following 
discussion, Ar = 0. The $ function for the ladder systems investigated in this paper 
is shown in figure [2l corresponding to a screened Frohlich interaction with Rg^ = 1- 
We also define the dimensionless interaction parameter A = Ep/W where W is the 
magnitude of the energy of the tight-binding electron (normally the half band-width 
zt). As we can see from equation U\ the energy shift when there is only one particle 
(polaron shift for i = i') is Ep = f^j = Thus A = ^wtfJ' ■ Substituting 

that definition into the atomic Hamiltonian, one obtains: 

a' ^ ' j 

The reason for introducting the new functions $ can immediately be seen, since they 
appear in the Hamiltonian as a ratio. Thus, in combination with A they give a universal 
definition of coupling in models with long-range hopping. 

I Inspection of equations [3] and [3] shows that electron number operators are unchanged on 
transformation, so the Coulomb part of the Hamiltonian is also unchanged 
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Figure 2. $ functions for (a) the staggered ladder and (b) the rectangular ladder 
(subscripts 1 and 2 correspond to the leg of the chain). Note that since an index is 
assigned to each unit cell, there is an offset of 1 index in the interaction funcion between 
paths 1 and 2, relative to the interaction function between paths 2 and 1 when the 
staggered ladder is simulated. 



2.2. Transforming the electron hopping term 

Substitution of equations ([3lll]) transforms the tight binding Hamiltonian in the following 
way: 

Htb = tii>c\ci> ^ f^i^ = ^ tii,c\ci> = ^ CTii^clci' (10) 
ii' ii' ii' 

where 

aw = tie exp gij{d] - dj)^ exp ^- ^ gi'jid] - dj)^ (11) 

i.e. we can see that the electron-phonon interaction in the transformed Hamiltonian is 
part of the hopping process, and that to some extent, one may regard the operators c 
as creating polarons (i.e. electrons and a phonon cloud at the same time). 
We now apply the following identity, 

(which is valid if C = [A, B] commutes with both A and B) with e'^e"^ = 1, which 
is always valid. Therefore, e^^e^ = e^+^ei^'^lf^^ in equation [TT] we may choose 
^ = Ej 9ijid] - dj) and 5 = - 9i'jid] - dj). So, A + B = Y.ji.9ij - 9i'j)id] - dj) 
and [A, B] = Since the commutator is a number, equation [T2] may be applied. 
The hopping operator becomes: 



an' = tii> exp ^- Y^igij - gi'j)dj + Y^{gij - gi'j)d]j 



(13) 



J 3 

Now, the identity may be used again, making the grouping, A = — J2ji9ij — 9i'j)dj 
and B = Y.ji.9ij - 9i'j)d\- Thus [A, B] = - Y.ji.9ij - 9i'j)i9ij - 9i'j), and 



"r-'"J 

□ E,- -9i' , )d] p-T.i i9i3 -Oil Mj 



o'ii' = tii' exp 



WX $0(2,2') 



uj \ $0(0,0) 
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This form is particularly useful, due to the order of the creation and annihilation 
operators. Thus, when one carries out the perturbation theory, calculation is reduced 



It's interesting to note that, when computed as an average over the atomic wavefunction, 
the hopping integral may be regarded as being modified by the EPI as tw —>■ 
tii' exp {—W X'jii' / uj) = iii', where •yui = 1 — $o(^! ^')/^o(0, 0). Such an approximation 
is valid in the anti-adiabatic limit, and will be revisited later in this paper. The band- 
narrowing factor was originally introduced by Tyablikov using an equations of motion 
scheme |24j . 

3. Quantum Monte-Carlo simulation 

The CTQMC algorithm presented here is an extension of a similar path-integral method 
for simulating the polaron problem [19] . An integration over phonon degrees of freedom 
following Feynman leads to an effective action, which is a functional of two polaron 
paths in imaginary time which form the bipolaron and is given by the following double 
integral when huP ^ 1 |22j , 



where the vector Ar = r(/5) — r(0) is the difference between the end points of one of 
the paths in the non-exchanged configuration (here u = hu/t and (3 = t/ksT). The 
indices i = 1,2 and j = 1,2 represent the fermion paths. \^(ri,r2) is an instantaneous 
Coulomb repulsion. The part of the action depending on Ar arises because the entire 
phonon subsystem at r = /? must also be shifted when there is a shift in the electron 
sub-system between the start and end configurations. The definition of Ar and other 
nomenclature for the CTQMC simulation of ladder systems are shown in figure [31 

From this starting point, the bipolaron is simulated using the Metropolis Monte- 
Carlo (MC) method. The electron paths are continuous in time with hopping events (or 
kinks) introduced or removed from the path with each MC step. Analytic integration is 
performed over sections of parallel paths. The ends of the two paths at r = and t = (3 
are related by an arbitrary translation, Ar. In contrast to the one-particle case, the 
fixing of the end configurations limits the update procedure to inserting and removing 
pairs of kinks and antikinks. We constrain particles to opposite legs of the ladder, which 





X 



5^($Ar[r,(r),r,(r')]-$o[r.(r),r,(r')]) 




(15) 
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Figure 3. Example path configuration on the ladder system, showing the notation 
used. Paths are separated by a vector b, and sites in the chain by a. The path 
configuration at r = is identical to that at r = /3 up to a shift Ar (i.e. the ends of 
the beginning and end of the paths are separated by A). Each path has a number of 
kinks Na/b+ and antikinks Na/b- 



corresponds to two species of charged particles. In such a system, there is no exchange 
between particles. Exchange and singlet triplet splitting from Quantum Monte-Carlo 
simulations is briefly discussed in this section, with an analytical discussion in section 
5. A full discussion of the QMC procedure for exchange is left to a future article. 

3.1. Binary updates 

In the path integral QMC with two paths, it is necessary to make two kink 
operations simulataneously to ensure that the end configurations remain identical 
up to a translation. There exist two classes of update. First, kink/antikink pair 
addition/removal on a single path is useful, since it always maintains the same end 
configurations up to a change in the interpath displacement. Second, kink pair additions 
on different paths are needed in the analysis of the bipolaron mass, where the t — (3 end 
configuration must be exactly equal to the r = end configuration up to a parallel shift. 
Within both subclasses, there arc two specific update types that satisfy the imaginary- 
time boundary conditions, as follows: 

(I) Two kinks of the same type 1 are added to or removed from two different paths. 
(II) A kink-antikink pair is added to or removed from one of the two paths. An antikink 
to kink 1 is a kink with the opposite direction —1. 

(III) A kink of type 1 is inserted into one path, and another kink of the same type 1 is 
removed from the same path (kink shift). 

(IV) A kink of type 1 is added to one path, and an antikink —1 is removed from the other 
path. 

An important property of the bipolaron system is that the type and time of added 
or removed kinks still does not define the new path unambiguously. Indeed, imagine a 
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kink of type 1 being inserted on a single path at time Ting. This could change the path 
in two different ways: Either the path at times r > Tins is shifted in the direction 1, or 
the path at times r < Tins is shifted in the antidirection 1. We refer to the former change 
as a top shift, and the latter as bottom shift. For the single-path (polaron) problem 
the distinction between the top and bottom shifts was not important because they are 
identical up to a translation of the entire path. This argument does not hold for two 
paths, since the resulting interpath distance in a binary update changes with shift type, 
thus a choice of shifts is an important part of the Monte Carlo update process. We 
proceed to derive the update probabihties for the Monte-Carlo scheme set out above. 

There is considerable flexibility in choosing the probabilities for adding and 
removing kinks. We choose an equal weighting scheme for choosing kinks, shifts and 
paths as follows: 

Ai Choose a kink type 1 from the list of all possible kinks, with equal probability 
Pi = l/Nk, where is the total number of kink types. Pi cancels on both sides of 
all balance equations considered below. 

Aii Anti-kink types are always determined from the kink type. 

Aiii Shift type (top or bottom) is chosen with equal probability P, = 1/2 independently 
for the two kinks. Pg also cancels on both sides of the balance equation. 

Aiv Assign path A with equal probability 1/2 from the two available paths. 

Av Assign path B as the other path. 

We will also choose kinks to add and remove according to the following weightings, 
although there arc some specific rules in the updates below to deal with cases where 
there are no kinks or antikinks of the chosen type on a path. 

Bi The probability density for kink time selection when adding a kink is always 

p{t)^1/(3, (0<t</3). 
Bii The probability for removing a kink of type 1 from path A in a configuration C is 
i/NAi{C) where Nai{C) is the number of kinks of type 1 on the path A. 

We note that the following is not the only set of possible rules, however, we consider 
this to be the most transparent method of choosing kinks to insert and remove. 

Update type (I): Addition/removal of two kinks to or from different paths Consider 
two configurations with two paths, C and D, where configuration D has two more 1 
kinks than C, one on the first path at time ri, and one on the second path at time T2. 
The balance equation is 

WiC) ■ Qa{C) ■ P{C ^D) = W{D) ■ Qr{D) ■ P{D -> C) . (16) 

The relative weight of configurations C and D is W{D)/W{C) = (tiAr)2e^(-^)-^(^). 
In order to approach the limit of continuous time, we rewrite the probability Qa of 
selecting two kinks at Ti and T2 to add to different paths given a configuration C as 
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QaO^, Ti, A; 1, T2, B\C) = QaO-, ti, A; 1, T2, i?|C)(Ar)^, where qa is the corresponding two- 
dimensional probabihty density. The probabihty of selecting the two specific kinks Ti 
and T2 to remove from the resultant configuration (D) to get back to C is a finite 
number Q^{\,Ti, A;1,T2, B\D), which is the probability of removing those two kinks 
given configuration D. Substituting these results into the balance equation, one may 
cancel the (Ar)^. Applying Metropolis, update probabilities are obtained as follows: 

P(C - i^) = min 1 1 ; ^i^Ra -4, 1, r,, ^^(^)_^(^) \ 
^ ^ I qR{l,n,A;l,T2,B\C) } ^ ^ 

P{D ^ C) = min |l ; gA(l, r,, A 1, r,, i^jC) | _ 

Thus we have obtained update probabilities that do not depend on the time 
discretisation, and we can immediately take the limit of continuous time. A similar 
approach can be taken for any of the update types I to IV. We now demonstrate all 
steps in the derivation of the first update probability as an example. 
The rules and resulting probabilities are as follows: 

(i) Choose kink types, shifts and paths according to rules Ai-Av. 

(ii) If the initial configuration has at least one 1 kink on path A and on path B, then 
removal of a pair is proposed with probability Pr = 1/2, and addition of a pair 
is proposed with probability Pa — 1/2. Otherwise, only pair addition can be 
attempted and Pa — 1- 

(iii) If pair addition is selected, times are selected to insert one kink on path A, and 
another on path B with independent equal probability density (rule Bi). 

(iv) If pair removal is chosen, then one candidate kink is selected with independent 
equal probability from each of paths A and B in configuration D (rule Bii) . 

Implementing these choices, one obtains qA{l,Ti, A;\,T2, B\C) — PiPg Pa{C) / P'^ from 
the combination of rules (i), (ii), (iii). Likewise, the combination of rules (i), (ii) 
and (iv) specifies that QR{\,ri,A;l,T2,B\D) = PiP^PR{D)/NiAiD)NiBiD). Rule (iii) 
leads to Pa{C) = 1/2 if Nai{C) > 1 and Nbi{C) > 1 and Pa{C) = 1 otherwise. 
Since configuration D always has sufficent kinks to make a removal, we always have 
Pr{D) = 1/2 

leading to the following acceptance rules: 



Pij(D)(ti/3)2e^(^)-^(^) 

Please note which configuration the kink numbers apply to. In the case of kink 
addition, initial configuration is C and final configuration is D. In the case 
of kink removal, the initial configuration is D and the final configuration is 
C. 
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Update type (II): Addition/removal of a kink-antikink pair to one path The general 
properties of this update type are similar to update type (I): the addition of a 
pair at times ri and T2 is characterized by a two-dimensional probability density 
qAii,Ti, A; —l,T2, A\C) while removal of a pair is charaterized by a finite number 
Qr(I, Ti, A; I, T2, A\D). Eqs. (fT7j) -( fT8l) still apply with in place of tf. Since ti = 
the acceptance rules are idential, and only Q and q differ. 
We use the following rules. 

(i) Choose kink types, shifts and paths according to rules Ai-Av. 

(ii) If the initial configuration has at least one 1 kink and one antikink {Nj^i > 1, 
Na-\ > 1), then removal of a pair is proposed with probability Pr = 1/2, and 
addition of a pair is proposed with probability Pa = 1/2. Otherwise, addition of a 
pair is proposed with probability P4 = 1- 

(iii) If pair addition is selected, kink and antikink insertion times are selected with 
independent equal probability density for insertion on path A. (Rule Bi) 

(iv) If pair removal is chosen, then one candidate kink and one candidate antikink are 
selected with independent equal probability from path A in configuration D (Rule 
Bii). 

One computes the update probabilities as before: 

f PR(D)(t^3)^e^^^^-^^^^^ 
P(addition) = P{C ^D) = inm\l; „ /J^,,, , \ (21) 



Pa{C)Na,{D)Na-i{D) 

Update type (III): Addition and removal of a kink to one path (kink shift) This 
update type does not change the number of kinks, and hence does not change the 
kinetic energy of the system. We define a configuration C, and a configuration D which 
is identical to C except that one kink has been shifted. To get from configuration C 
to D, a kink is removed from path A at time ti and is reinserted in the path at time 
T2. Since C and D have equal total kink number, the ratio of statistical weights is 
W{D)/W{C) = e^(^)-^(C'). From detailed balance and Metropolis: 



P(C ^ - min 1 1 ■ Qr(1> ^2, A\D)qAil, n, A\D) a(d)-A(c) ] 
^ ^ \ ' QK{\,r,,A\C)qp.{\,r2,A\C)' ] 



(23) 



There is only one update rule, since we can get from D and C using exactly the same 
process as going from C to D. All attributes of the kinetic energy have dropped out 
from the equations. The acceptance rules are determined solely by the electron-phonon 
interaction, as expected for this update type. 

In many practical situations it is reasonable to choose the functions q and Q 
independent of time r. In this case the above expressions simplify significantly. In 
particular, consider the following set of update rules. 

(i) Choose kink types, shifts and paths according to the general rules Ai-Av. 
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(ii) If the path has no 1 kinks, Ni — 0, the update attempt is aborted. Otherwise, 
A^i > 1, so propose an 1 kink for removal with equal probability l/Ni{C) (rule Bii). 

(iii) Propose time for new 1 kink with constant probability density p{t) — 1/(3 (rule Bi). 

These rules result in cancelation of Qremove and gfadd from the above equations, which 
reduce to 

P(C ^ L>) = min {1 ; e^(^)-^(^)} , (24) 

Update type (IV): Kink addition to one path and antikink removal from the other path 
In this update, a kink 1 is added to one path and an antikink —1 is removed from the 
other path. As a result, the (3 end configuration shifts as a whole by I. In a reciprocal 
process, a kink 1 is removed from one path, an antikink —1 is inserted in the other path, 
and the j3 configuration shifts by —1. 

The ratio of weights is W{D)/W{C) = {t-i/ti)e^^^^-^^'^\ The balance equation is 
satisfied by the following solution 

P(C^D)-mmll- -4|/^)gA(-l; r^, B\D) \ .35) 

F(C^D)-mm|l, Q^^_y.^^.B\c)<i^^y^n,A\cf ] 

Since we can obtain the inverse process by changing kink 1 for its antikink, the update 
probability P{D — > C) is not necessary since we always choose 1 from all kink types. 
As the simplest implementation, the following rules are used, 

(i) Choose kink types, shifts and paths according to rules Ai-Av. 

(ii) If path B has no antikinks —1, Nb-\ = 0, then the update attempt is aborted. 
Otherwise, Nb-i > 1, so an antikink is proposed for removal from path B with 
equal probability (rule Bii). 

(iii) The time location for kink insertion on path A is proposed with equal probability 
density (rule Bi). 

With these assumptions, the acceptance probability takes the form: 

P{C ^ D) = min { 1 ; ^^^^e^^^^"^^^^ } , (26) 

3.2. Estimators 

When our Monte-Carlo scheme has reached equilibrium, we make a series of 
measurements of physical properties. The ground state energy is: 



eo = — lim 

/?— >oo 



(27) 



where Ng is the number of kinks of type s, and angular brackets denote ensemble 
averaging. The number of phonons is given by: 

1 / dA \ 
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where the derivative is taken keeping Xu constant. The polaron band energy spectrum 
can be computed from: 



Ck — Co = — hm — ln(cos(k • Ar)) , 

/3^oo fj 



(29) 



where k is the quasi momentum. By expanding this expression in small k, the i-th 
component of the inverse effective mass is obtained as 



hm -5^((Ar,)2) 



(30) 



Thus the inverse effective mass is the diffusion coefficient of the polaron path in the limit 
of the infinitely long "diffusion time" /3. The bipolaron radius is the average distance 
between paths, 



Rbp — 




Ari2(T)2dT 



(31) 



Finally, the mass isotope coefficient, a 

CO 1 



d\nm*/dlnM, is calculated as follows 



am* = lim ,, , 

2 ((Ari)2) 



- ((Ar,)^) ( ^ 



/dA 



(32) 



3.3. General Monte-Carlo considerations 



There are certain aspects of good practice for quantum Monte-Carlo simulations that 
we adhere to here. As always for Montc-Carlo simulations a random number generator 
with sufficient period is used. Measurements arc performed every few steps to avoid 
unnecessary correlations in results (the aim here is to spend no longer measuring than 
simulating, since time correlated results do not make a large contribution to more 
accurate measurement). Careful blocking analysis with large blocking sizes Nb is 
performed to determine accurate error bars. To avoid anomalous error bars caused 
by long time correlations, we compare error bars computed with two block sizes Nb and 
2Nb. 



3.4- Exchanges 

Exchanges are significantly more complicated, with several possibilities for the exchange 
update involving inserting and removing kinks. In the exchanged configuration, there 
are an additional 4 update rules, and there is also an ambiguous configuration where 
both paths have the same start and end points, which leads to some small additional 
modifications. We defer a full discussion of exchange update rules to a later paper. On 
our ladder models, exchanges are not required, since electrons sit on opposite legs. 



3.5. Singlet-triplet splitting in the Monte-Carlo method 

A consequence of exchange is that singlet and triplet states arc not degenerate. We can 
see the singlet-triplet mass difference as a consequence of interference between paths in 
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the Monte Carlo simulation. We take a simplified one- dimensional example to illustrate 
the mechanism. Consider a bipolaron of separation R propagating from the sites {0, R} 
at r = to { Ar, R + Ar} at t = (3. We assume that the weight w of a single-electron 
path is a monotonically decreasing function of the number n of kinks in the path, and 
the paths with the smallest number of kinks dominate. This is likely to be valid in the 
strong-coupling limit. We can therefore write the weight of as path as a monotonically 
rapidly decreasing function w{d), where d is the distance between endpoints. We also 
neglect the interaction between paths. 

Consider first the case of periodic boundary conditions Ar = 0. The weights of the 
singlet and triplet bipolaron paths are respectively the sum and the difference of the 
direct and exchange paths. 

w,(0) ^ w{0)w{0) + w{R)w{R), (33) 
WtiO) ^ w{0)w{0) - w{R)w{R). (34) 

Here if w{R) ^ w^O), the singlet and triplet weights are dominated by the direct path 
and are nearly equal. 

Now consider a twist larger than the bipolaron radius, Ar > i? > 0. The singlet 
and triplet weights are dominated by the shortest paths: 

w,{0) = w{Ar)w{Ar) + w{Ar + R)w{Ar - R), (35) 
wt{0) = w{Ar)w{Ar) - w{Ar + R)w{Ar - R). (36) 

In this case the total number of kinks in either the direct or the exchange path is 
the same (2Ar), so there will be substantial cancellation in the triplet case. Thus the 
diffusion coefficient of the triplet bipolaron will be smaller, and the effective mass larger, 
than that of the singlet. 

4. Polarons on triangular and rectangular systems 

We briefly discuss the lattice dependent features of the polaron problem here. For more 
detailed discussion of the polaron problem, the reader is directed towards our papers 
on this subject [21], Eni [19], and to the paper by Kornilovitch in this issue [26]. In this 
section, we specifically discuss the changes to the masses of polaron moving on trangular 
and square lattices as the screening length Rsc is varied. The polaron mass forms part 
of the argument later in this article. 

Screened Frohlich electron-phonon interactions were simulated by Spencer et al. 
[20] , who demonstrated a continuous crossover between the Frohlich and Holstein limits 
on the chain. In particular, the mass of the particle is found to be light down to quite 
small screening radii, consistent with results by Bonca and Trugman [18] for nearest- 
neighbour electron-phonon interactions. 

We have previously computed the properties of polarons on several Bravais lattices, 
showing that the effects of the lattice type on the properties of the polaron are 'washed 
out' by long range interactions [21]. Figure H] shows the effective mass of the dicrete 
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Figure 4. Inverse effective mass of tlie discrete Frohlich polaron on square and 
triangular lattices, u/t = 1, and A is varied. Triangular lattices are shown on the 
left, and square lattices are shown on the right. From top to bottom, the screening 
radius of the interaction is decreased, with the top graphs showing Frohlich interaction 
Rsc oo, middle graphs screened interaction R^c — 1 and the bottom graphs the 
Holstein interaction Rsc — 0. Frohlich polarons are significantly lighter than their 
Holstein conterparts because the long range interaction leads to pre-distortions of the 
lattice before a hopping. 



Frohlich polaron on square and triangular lattices. Frohlich polarons are significantly 
lighter than their Holstein counterparts, due to the long range interaction. We have 
also shown that the overriding factor for the properties of Holstein polarons is the 
number of nearest neighbours in the lattice, and not the dimensionality [21]. Since we 
are interested in long-range interactions in this paper, then it suffices to note that the 
masses of Frohlich polarons are of the same order of magnitude on all lattice types, and 
that they are extremely light [21]. The properties of realistic screened interactions lie 
somewhere inbetween, but such polarons remain light down to quite small interaction 
ranges of the order of a lattice spacing [20] . 
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5. Bipolarons on ladder systems 

The properties of bipolarons are significantly more complicated than those of polarons. 
One expects bound Holstein bipolarons in strongly correlated materials to be extremely 
heavy, since a very large attractive potential is needed to overcome the repulsion. This 
is not true for longer range interactions. According to work by Bonca and Trugman |18j . 
even bipolarons with nearest neighbour interactions are significantly more mobile than 
their Holstein counterparts (indeed, one normally expects the mass of a bipolaron to 
scale like the squared mass of the polaron, so polarons which are an order of magnitude 
lighter than Holstein ones, may become bipolarons which are two orders of magnitude 
lighter). Another extremely interesting proposition is the role of geometry on the 
bipolaron mass. When bipolarons are bound on nearest neighbour sites, and another 
degenerate state may be reached in a single hopping event, the leading correction to 
the atomic Hamiltonian is first order in the hopping term, and not second order as one 
might expect for the Holstein model, which leads to a bipolaron with a mass that is of 
similar order of magnitude to the polaron mass (i.e. a superlight bipolaron) [9j. Such 
systems may be realised on triangular lattices, or on lattices with large next-nearest 
neighbour hopping. We have recently extended our Quantum Monte-Carlo algorithm 
to explore this type of bipolaron, leading to similar conclusions [22j . In this paper, we 
discuss some of these extensions to the algorithm to look at two types of ladder system 
shown in figure [1] 



5.1. Weak and strong coupling 

Since the particles on the different legs of the ladder cannot exchange the very weak- 
couping limit is not well bound, consisting of two large polarons. As such, the weak- 
coupling perturbation theory can be made about the unbound state. As will become 
apparent when we show the quantum Monte-Carlo results, the number of phonons 
is small for the weakly bound states, especially in the anti-adiabatic limit. The 
perturbation theory in the electron-phonon coupling term only excites single phonons, 
so if the number of phonons becomes too large, the theory fails. In the presence of strong 
on-site repulsion, the bipolaron is not bound at zero coupling. In general, if there is 
a bound state for A — 0^, then this perturbation expansion fails. The expansion is 
written as follows: 

_ (0) o A^iy 1 ^ l/qP 

^"--^'^ '^%io:o)N^w{k:^)^ ^^^^ 

W{k,ci) = e^^l^ + uj-e^^\ (38) 
/q = $^/m(0)e-'^-, (39) 

m 

where N is the number of momentum states. Thus the ground state polaron energy (at 
k = 0) is eQ^'' = —W + XWT^g{uj), which defines a dimensionless coefficient Teq- 
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There is no general analytic solution for the second order perturbation theory, but 
values may be computed using numerical integration. The number of phonons, isotope 
exponent and inverse mass may also be written as, Nph = 2\Tn{}tio), a = 2\Ta{fkj) 
and 2mo/m** = 1 — rm(^)A. The weak coupling limits for the polaron are discussed 
in references [20l |2T] amoung others. A is defined as before. 

Aspects of the strong coupling limit are easy to compute from the path integral 
formalism. In the very strong coupling limit, the most common configurations are 
straight paths, since the action is proportional to A and that configuration gives the 
smallest possible action. Thus the strong coupling action reads: 

= ^0(0,0)^ ^ E^o(.j) (40) 

For a polaron, the total energy for a self-interacting strong-coupling system (i.e. a 
straight path) is E = —WX, where W is the half band width, — ek=o- Computing 
the action for the straight path configuration, one obtains the energy as -^strong = 
-Mstrong5/?|/3^oo- Leading to, 

^'t?ong = -2W^A[1 + $o(0, b)/$o(0, 0)] (41) 

Where b is a nearest neighbour vector between chains indicating the point of closest 
approach. This result is of little suprise, since there are 2 polarons in a bipolaron, with 
E = —WX each, and 2 inter-polaron binding terms with E = —WX^{0, b)/$(0, 0). We 
may also compute the number of phonons associated with the bipolaron in a similar 
manner using equation EH] with ^strong, which we obtain as: 

^^s'^tLg = ^WX[l + <f (0, b)/$(0, 0)]/^ (42) 

As we have discussed, the $ functions for the ladder systems with Rsc = 1 are 
plotted in fig. d Numerical values for $(0,b), $(0,0) and 1 + $(0, b)/$(0, 0) are 
shown in table [H from which numerical values for the strong coupling behaviour are 
computed. We discuss the region of validity of these results when we compute the exact 
QMC results. In particular, if we plot appropriate ratios, such as the ratio Nphu/X, we 
expect saturation at strong coupling, consistent with table [H 

Typically, it has been discussed that an expansion of the Lang-Firsov transformed 
Hamiltonian in the polaron hopping can be used at strong coupling. The perturbation 
expansion to second order may be written as, 

Etot = E^, + mH^m)^, + V (43) 

J ■' 

The values for the strong coupling energy may be computed from the leading term 
of this expansion, and have the same form as the energy computed from the straight 
paths. This is not suprising, since the corrections to the energy are at least as small as 
t. 

There is a subtle point associated with phonon numbers. By examination, the 
ground state of the phonon subsystem in the transformed Hamiltonian is d}jdj = 
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Lattice 


$(0,0) 


$(0,b) 


1 , <S>(0,b} 
^ $(0,0) 


E/\ 


ujNph/X 


Staggered 


1.06896 


0.30034 


1.28096 


-5.12384 


5.12384 


Rectangular 


1.05554 


0.284832 


1.26984 


-5.07936 


5.07936 



Table 1. Strong coupling behaviour of the staggered and rectangular lattices. 



rij = 0. In order to determine the total number of phonons in the true ground 
state of the atomic Hamiltonian, one must transform back to the regular wavefunction. 
One may either transform the wavefunction, or the phonon part of the Hamiltonian. 
Transforming the phonon part of the Hamiltonian is easy, and one obtains, 

\Hph\(j>LF) = {(p\Hph\<f>) = ujNph (44) 



thus 



Nph = {(I>lf\ ^{d] + y^^gijni){dj + '^gi:jni>)\(pLF)/uJ (45) 

j i V 

= 9ij9v jniUj' I / (46) 

ii'j 

the expectation value of the phonons occupation may be rewritten in terms of A and 
thus: 

Npu = -— (47) 

OJ 

(Note the minus sign in front of the equation - phonon number is positive). 

Thus one determines that the total number of phonons in the bipolaron case is: 



00 V ^o(0,0) 
(N.B. A similar argument can be used for polarons.) 

Unfortunately, the computation of the mass is not so simple. The mass is very 
sensitive to the exact form of the renormalised hopping, and the leading order of the 
perturbation expansion varies with lattice type. For the superlight small bipolarons 
on triangular plaqucttcs discussed here, there is a leading term with order t, and 
for rectangular systems, the leading order is t^. Examination of the path integrals 
demonstrates that the mass is computed from a parallel shift. The first perturbation of 
the paths from the straight line in the atomic limit is the insertion of two parallel kinks 
(one on each path). In order for the mass to be genuinely first order in the polaron 
hopping, we would require updates with only one kink. Thus, there are significant 
contributions from the second order term, including a cancellation between the orders 
as A — > CX3. In fact, the expansion of the Lang-Firsov transformed Hamiltonian to 1st 
order in the hopping turns out to correspond to the extreme anti-adiabatic limit, which 
we now discuss. 



5.2. Anti-adiabatic approximation 

The bipolaron dispersion can be evaluated analytically in the anti-adiabatic strong- 
coupling limit {huj 3> t, A 3> 1), in which case the problem can be reduced to that of 
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a rigid dimer. In this limit the Lang-Firsov transformation [25] ehminates the phonons 
from the Coulomb-Frohhch Hamiltonian ([1]). In combination with an averaging over 
phonons consistent with the anti-adiabatic approximation (if the phonon frequency is 
very large, there are no real phonons) one obtains the following effective UV Hamiltonian 

H — — ^ ^ ^nn'Cno-'^n'o" "~ ^ ] cJ^^Cno- 
nn'cr ncr 

+ ^ J] ^nT^nTCniCni + ' 5Z ^nn'cLcn<xC|,,^,Cn/^/ . (49) 

n nn' aa' 

where the primed sum excludes self-interaction, and the renormalised on-site and inter- 
site interactions are 

ij^m^-WX (50) 

and 

- —2 $o(0,0) ^^^^ 

respectively. (There is no retardation in the anti-adiabatic limit.) Now suppose the 
on-site repulsion U ^ t is large and repulsive and the inter-site interaction has a well- 
defined minimum Vmin < at some separation. Let A^n be the set of sites of n' at that 
separation from n: 

iVn = {n' : = KniJ (52) 

and 

Kn' - Klin > i for n ^ iVn. (53) 

We shall call A^n the "neighbours" of n. In general, A^n need not be the hopping 
neighbours {n' : tnn' 7^ 0}. The low-energy sector for two electrons corresponds to 
dimer states (bipolarons) in which the electrons are on neighbouring sites; the energies 
in this sector are Vmm + 0(t). The energy cost of internal excitations of the bipolarons 
introduces a gap. 

We can now sharply distinguish two types of bipolaron motion: "crab-like", in 
which the constituent polarons remain neighbours, and "crawler" , which requires virtual 
transitions out of the low-energy sector. The crab-like bipolaron bandwidth will be 
0(t), while the crawler contributions will be 0{t^). We shall for the present purposes 
project onto the low energy sector and hence ignore the resulting higher-order terms in 
t, immobilising crawlers. For simplicity we drop the tildes from the notation and absorb 
the polaron shift —Ep into the chemical potential. 

If a lattice A has L sites with a mean number v of neighbours per site, then the 
single-polaron Hilbert space is 2L-dimensional, the two-polaron Hilbert space is 4L^- 
dimensional and the crab bipolaron Hilbert space is just 4i/L-dimensional. We can 
further reduce to one singlet and three triplet spaces, each of dimensionality vL. The 
singlet bipolaron space is 



S = span (|n T n' i) + |n' t n i)) : n G A, n' G A^nj 



(54) 
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and the = 1, 0, —1 sectors of the triplet bipolaron space are 

Ti = span {|n t n' t) : n e A, n' G Nn} (55) 
To =span|^(|ntn'i)-|n'tni)) :nGA,n'GiVnl (56) 



V2 

T_i = span {|n i n' i) : n G A, n' G N^} . (57) 

This enables us to write the low-energy effective Hamiltonian of the dimers in each 
sector as a tight-binding Hamiltonian on the dimer lattice constructed in the following 
way: a node is placed on the line joining neighbours in the lattice A. If n' and n" 
are both hopping neighbours of n (tn'n" 7^ 0), then the dimer can hop from nn' to 
nn". A bond is then drawn between the two nodes on the dimer lattice with hopping 
integral tn'n" in the singlet sector and — tn'n" in the triplet sector. This sign change 
ensures the correct exchange symmetry for closed paths on odd-membered rings: as a 
dimer completes one cycle of an odd-membered ring its end-points are interchanged. 
This can lead to a dramatic difference between singlet and triplet bipolaron masses on 
a non-bipartite lattice, as we shall see. 

5.3. Ladders 

In the staggered ladder depicted in Fig [1] the neighbours of a site on one chain are the 
two adjacent sites on the opposite chain, while the hopping neighbours of a site are 
along the same chain. The corresponding dimer lattice is a one- dimensional chain with 
hopping ±t and two sites per unit cell. As there are no exchange paths, the sign of the 
hopping can be gauged away and there is no singlet-triplet splitting. The polaron and 
bipolaron dispersions are therefore respectively 

Epoi{k) = — i cos ka, (58) 

Eupik)= ±i cos— . (59) 

The bipolaron effective mass is therefore four times that of the polaron (as previously 
reported [9]): 

_2 d'^Epoiik) 



m** = 



dk^ 
d'E^pik) 



(60) 



k=0 



dk^ 

This result is remarkable if one considers the standard strong coupling result for the 
mass of the bipolaron. In the rectangular ladder the requirement for a virtual internal 
excitation of the bipolaron in the crawler dynamics would lead to a hopping ©(t^) and 
hence mass 0{{m*Y). The staggered ladder with long range electron phonon attraction 
has two degenerate nearest neighbour bound states (as summarised in figure 12), so 
no intermediate state is required when the particle hops. It is clearly important that 
the electrons are bound one lattice spacing apart to take advantage of this effect, but 
this can easily be achieved in the presence of a strong site-local Coulomb repulsion 
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(so that the energy is not at a minimum when both particles are on the same site). 
This spaced minimum is clearly extremely important when real lattices beyond the 
toy ladder models presented here are considered. Some details of real lattices will be 
discussed later in this paper, but to briefly summarise, in order to obtain this special 
kind of bipolaron, (a) The lattice must have several degenerate near-neighbour bound 
states which can be transformed from one degenerate state to another via a single hop 
(b) Strong Coulomb repulsion is required to stop a unique single-site bound state from 
forming between the polarons, and (c) long range attraction is required so that the 
minimum in the potential function (attraction + repulsion) is at approximately one 
lattice spacing. This information is summarised in figures [12] and [13] in the conclusion 
to this paper. 

5.4. QMC 

We compute QMC results for bipolarons moving on staggered and rectangular ladders 
with period 1000 for a range of A and u, including the total energy, inverse bipolaron 
mass, bipolaron radius, mass isotope exponent and the number of phonons in the 
bipolaron cloud. Figure [5] shows the total energy of the bipolaron for (a) the rectangular 
ladder and (b) the staggered ladder. A slight change in gradient between weak and 
strong coupling is just discernable, demonstrating (as we shall see in the coming figures) 
that the staggered ladder reaches the strong coupling limit at significantly lower A than 
the rectangular ladder. Strong coupling results from the previous section agree well. We 
see that there are no significant differences between the total energy of the bipolarons 
on the staggered and rectangular ladders, although the strong coupling limit is reached 
for slightly lower A in the case of the staggered ladder. 

If one is to reach a superconducting state via the Bose-Einstein condensation of 
bipolarons, there are two conditions. First, the bipolaron pair must be light, and second, 
the pairing radius must be small. We demonstrate the differences between the inverse 
masses on the staggered and rectangular ladders in figures [6](a) and (b), which show the 
inverse mass of the bipolaron for a number of different A and u/t. There is more than 
an order of magnitude difference between the mass of the bipolaron on the staggered 
ladder, and that on the rectangular ladder over significant regions of the parameter 
space. In fact, the magnitude of the bipolaron mass turns out to be of similar size to 
that of the polaron mass over a wide range of the parameter space. The mass is inversely 
proportional to the transition temperature of the BEC, so a small mass is essential to 
obtain a decent Tc- 

We have demonstrated that one of the precursors for a Bose-Einstein condensate of 
pairs above the mK range may be met on the staggered ladder arrangement, but we also 
require small pairs, with non-overlapping wavefunctions. In figure [7] we show how the 
the size of the bipolaron varies as A and u/t are varied. Not only is the bipolaron on the 
staggered ladder lighter than that on the rectangular ladder, it is also has a significantly 
smaller radius than the bipolaron on the rectangular ladder, making it a much better 
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Figure 5. Total energy of the bipolaron for (a) the rectangular ladder and (b) the 
staggered ladder. A slight change in gradient between weak and strong coupling is just 
discernable, demonstrating (as we shall see in the coming figures) that the staggered 
ladder reaches the strong coupling limit at significantly lower A than the rectangular 
ladder. 

prospect for Bose-Einstein condensation. 

Since the result that the bipolaron mass is proportional to the polaron mass via 
a numerical value relies on the anti-adiabatic limit, where the phonon frequency is 
very large, so that phonons may not be excited, we investigate the number of phonons 
asociated with the bipolaron in figure [HI The result is weighted by phonon frequency 
and electron phonon coupling, so that the strong coupling result can clearly be seen. 
Again, we can see that the strong-coupling limit is reached at significantly lower A on 
the staggered ladder than on the rectangular ladder. Strong coupling results from the 
previous section agree well, with the ratio urNph/ A approaching the numerical value given 
in table [H The anti-adiabatic limit can clearly be identified as regions on the graph 
where the phonon number approaches zero (i.e. the u > X quadrant of the parameter 
space), consistent with results for the mass. 

We show the mass isotope exponent of the bipolaron in figure O Again, we can 
see that the strong coupling limit is achieved at significantly lower A on the staggered 
ladder, compared with the rectangular ladder. The isotope exponent is also smaller on 
the staggered ladder, demonstrating a much smaller range of mass from weak to strong 
coupling. 

Finally, in figure [THl we show example path configurations on (a) the rectangular 
ladder and (b) the staggered ladder. Hopping events on different paths on the 
rectangular ladder are very closely correlated on the imaginary time axis. On the 
staggered ladder, there are two degenerate configurations, and paths are just as likely to 
sit on either of the two neighbouring sites, significantly reducing the correlation between 
kinks, and increasing the probability that kink pairs can be inserted. It is this that leads 
to significantly lighter bipolarons on the staggered ladder. 
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Figure 6. Inverse mass of the bipolaron for (a) the rectangular ladder and (b) the 
staggered ladder. There is more than an order of magnitude difference between the 
mass of the bipolaron on the staggered ladder, and that on the rectangular ladder over 
significant regions of the parameter space. Bipolaron masses on the staggered ladder 
have recently been shown to have a value commensurate with the polaron mass. 




Figure 7. Size of the bipolaron for (a) the rectangular ladder and (b) the staggered 
ladder. Not only is the bipolaron on the staggered ladder lighter than that on the 
rectangular ladder, its size is also smaller for equivalent A. 




Figure 8. Number of phonons associated with the bipolaron for (a) the rectangular 
ladder and (b) the staggered ladder, weighted by phonon frequency and electron 
phonon coupling. Again, we can see that the strong coupling limit is achieved at 
significantly lower A on the staggered ladder. It is interesting to note that the states 
with large omega have a smaller number of phonons. This is consistent with the anti- 
adiabatic approximation: Whc;n phonon frequency is large, creating phonons becomes 
difficult, and phonon wavefunction is the vaccum state. Then the phonon problem 
maps onto a UV model (discussed later). 
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Figure 9. Mass isotope exponent of tlie bipolaron for (a) tlie rectangular ladder and 
(b) the staggered ladder. Again, we can see that the strong coupling limit is achieved 
at significantly lower A on the staggered ladder. The isotope exponent is also smaller 
on the staggered ladder, demonstrating a much smaller range of mass from weak to 
strong coupling. 
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Figure 10. Example path configurations on (a) the rectangular ladder and (b) the 
staggered ladder. Hopping events on different paths on the rectangular ladder are very 
correlated in time. On the staggered ladder, there arc two degenerate configurations, 
and paths are just as likely to sit on either of the two neighbouring sites, significantly 
reducing the correlation between kinks, {oj/t = A = 4 in both cases.) 



6. Beyond ladders: Other superlight systems 

On the ladder systems, the electrons were held on neighbouring legs of the ladders. 
This is partially consistent with a very strong local Coulomb repulsion or Hubbard U. 
It is also possible to examine the effects of a very strong or even infinite Hubbard U 
on the masses of bipolarons bound via long range attraction on other low dimensional 
systems. As we have seen by examining the ladder systems, there are 3 requirements 
for superlight bipolarons (1) Electrons are not allowed to bind on a single site (2) There 
are at least two degenerate configurations of electrons sitting on neighbouring sites and 
(3) There are single hopping events which transform one configuration into another 
degenerate configurations. Then, hopping of the bipolaron is first order in the hopping 
of the polaron. Condition (2) is satisifed by long range electron-phonon intereraction, 
and there are several tight binding lattices which also satisfy conditions (2) and (3), and 
a strong site-local Coulomb repulsion satisfies condition (1). We discuss bipolarons on 
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these lattices in the anti-adiabatic approximation in this section. 



6.1. Triangular molecule 



The staggered ladder system discussed earlier in this paper can be considered to be 
made up of triangular plaquettes, so the first logical step to looking beyond the ladder 
systems is to analyse the physics of a single plaquette. If hopping is allowed between 
all the sites on the plaquette, then one may consider exchange effects, which lead to 
singlet and triplet bound states. To see exchange effects, wc require lattices with odd- 
membered rings. As the simplest example, consider three sites 0, 1,2, all of which are 
neighbours, with hopping t > 0. The polaron Hamiltonian is 

/ -i -i\ 



-t 

-i -i 
-i -i 



(62) 



The dimcr lattice is constructed by placing a node at the centre of each bond. The 
singlet bipolaron Hamiltonian is therefore 



y -i -i 
with eigenvalues {Kin — 2i, Kin + 
/ Kin i 



(63) 



K„„ ^ 

i, Kin + i}, and the triplet bipolaron Hamiltonian is 
i \ 



V 



K 



t 



t K 



(64) 



with eigenvalues {Kiin — i, Kin — i, Kin + 2f}. (An alternative choice of basis would 
be the six- dimensional unsymmetrised Sz = basis {|i f j J,)}. This would transform 
the triangle into a six- member ring. Such a basis becomes unwieldy for typical lattices). 
Clearly, the properties of the plaquette arc defined by single polaron hopping events, 
since there arc no eigenvalues with 0{t^) terms. It is this property which makes lattices 
which can be constructed from triangular plaquettes a good starting place when looking 
for small superlight bipolarons. 



6.2. Triangular lattice 

Let us consider a triangular lattice with nearest-neighbour hopping t. In the anti- 
adiabatic approximation, this is replaced by t. The polaron band structure (ignoring 
polaron shift) is -E'(k) = — 2fC(k), where we have defined 

C(k) = cos A;a;a + cos — \/3kya/2^ 

+ cos (k^a/2 + VSkya/2) , (65) 
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Figure 11. The triangular lattice (thin lines) with a dimer state (circle) at the mid- 
point of each bond. The polaron and bipolaron states are indicated on one triangle, 
with spin indices suppressed. The dimer lattice (thick lines) is a kagome lattice. 



a being the lattice parameter. Expanding near the F point gives 

E{k) = -6i+h^aH+0{k*) 
The polaron effective mass is therefore 



m 



?>ta? 



(66) 



(67) 



By placing a node on each bond in the lattice, we see that the resulting dimer lattice 
is a kagome lattice. The dispersion is easily found as follows (see also ref [25]) by 
diagonalising the secular matrix 

1 + 7* 1 + /? 
1 + 7 1 + a* 



H{\^) = V^i^ + t 







with 



a = exp (^—ikxa/2 + iVSkya/2^ 

(3 = exp (^—ikxa/2 — iVskya/2^ 
7 = exp {ikxO) . 



(68) 

(69) 

(70) 
(71) 



The sign in (1681) is —t for the singlet and +t for the triplet. There are three bipolaron 
bands with no gaps: 



^l(k) = Knin ± t (^-1 - V3 + 2C(k) 

^2(k) = Knin ± i (-1 + V3 + 2C(k) 
^3(k) = Knin±2t. 



(72) 

(73) 
(74) 
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with sign +t for the singlet and —t for the triplet. 
The lowest singlet band is 

^i(k) = K,in - 4t + Oik"), (75) 

with effective mass 

m** = ^ = 6m*. (76) 

We notice that this increases linearly with the polaron mass, indicating that crab-like 
bipolarons can be relatively light. 

The lowest triplet band is the flat band -E'3(k) = V — 2t. This implies that a triplet 
crab-like bipolaron has infinite mass on a triangular lattice. Once crawler motion is 
permitted, the effective mass is expected to be finite but proportional to (m*)^. 



6.3. Lattices with long range hopping 

Lattices with triangular components are not the only ones where the ability to 
move between degenerate paired states can be utilised to lead to small masses. For 
example, if we introduce next-nearest-neighbour hopping to a linear chain, then, in 
the anti-adiabatic approximation, we obtain a bipolaron mass, 2(2a)^t' and a polaron 
mass 2aH where i = texp{-WX{l - $(0, a)/$(0, 0))/u;) and i' = t'exp(-M^A(l - 
$(0, 2a)/$(0, 0))/ti;). While this does not lead to an exact cancellation of the exponents, 
the bipolaron to polaron mass ratio is linear in t'/t, and the bipolaron is therefore 
expected to be light. 

Another suggestion for bipolarons with light mass comes from the cuprate lattice. 
It has often been suggested that the tight-binding structure of the cuprates is a square 
lattice with nearest-neighbour hopping t and next-nearest-neighbour hopping t'. In the 
presence of very stong Coulomb repulsion U, one may determine that the ground state 
of the Lang-Firsov transformed Hamiltonian consists of nearest neighbour pairs with 
degeneracy four. If one applies the anti-adiabatic approximation to such a lattice, one 
determines that the hopping of such a bipolaron is first order in t'. We obtain a bipolaron 
mass, 2{y/2a)H' and a polaron mass 2a^t where i = texp{-WX{l - $(0, a)/<l>(0, 0))/^^;) 
and i' = t'exp{-WX{l - $(0, a - i?^°a)/<l>(0, 0))/cj). Here is the rotation operator. 
Again, while we don't get an exact cancellation of the exponents, we obtain a resulting 
mass proportional to t'/t, and light bipolarons are likely Such pairing is an 
appealing prospect, especially since the Hague has demonstrated the potential for d- 
wave superconductivity mediated by Holstein-like electron-phonon interactions in the 
intermediate coupling limit [29]. 



§ It is interesting to note that these lattices share a common feature with triangular lattices, namely 
that it is possible for a single particle to hop 3 times to return to the point of origin, indicating a more 
universal prospect for this phenomenon 
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Staggered ladder 




Figure 12. Pictorial demonstration of lattices which are expected to lead to 
superlight and light bipolarons when long range electron-phonon interaction is present 
in combination with a strong site-local Coulomb repulsion (Hubbard U) which stops 
on-site binding. All degenerate nearest-neighbour bound states are shown as the thick 
lines. The lightest bipolarons are expected on lattices with a trangular component, 
where hopping between degenerate states can be realised with the same single hop 
as the polaron (shown as arrows) . Thus the bipolaron moves in the first order of the 
polaron mass. We also show potential single hops for light bipolarons on the chain, and 
the square lattice, which can be realised if next-nearest-neighbour hopping is present. 
The superlight bipolarons are contrasted with traditional (Holstein) on-site bipolarons, 
which move with two hopping events, the first one breaking the bound state (step ii), 
and the second reforming it (step iii), thus such on-site (crawler) bipolarons cannot be 
simultaneously small (well bound) and mobile. 



7. Conclusions 

The bipolaronic extension [1] of the BCS theory towards the strong interaction between 
electrons and ion vibrations proved that the Cooper pairing in momentum space [Tj and 
the Ogg-Schafroth real-space pairing [271 EH] ci-re two extreme limits of the same problem. 
For a very strong electron-phonon coupling, polarons become self-trapped on a single 
lattice site and bipolarons are on-site singlets. In the Holstein model of the zero-range 
EPl their mass appears only in the second order of polaron hopping [18] , so that on-site 
bipolarons are very heavy. This estimate led some authors to the conclusion that the 
formation of itinerant small polarons and bipolarons in real materials is unlikely |30j . 
and high-temperature bipolaronic superconductivity is impossible [3T] . 

In fact, we have demonstrated here that small but light bipolarons could exist for 
realistic values of the finite-range EPI with high-frequency optical phonons in staggered 
ladder systems. Small light bipolarons are an essential precursor to high-temperature 
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Figure 13. Pictorial demonstration of the interpolaron effective interaction which can 
lead to superlight bipolarons (cross-section of anti-adiabatic limit of triangular lattice 
with A = 1). A strong on-site Coulomb repulsion (Hubbard U) stops on-site binding, 
leading to binding between polarons on neighbouring sites when there is a long range 
phonon mediated attraction. In this case, the potential well should lead to binding 
of small bipolarons on the order of one lattice site. The long range Frohlich electron- 
phonon attraction may be found on quasi-2D lattices, where ions oscillate above planes 
(as described in reference [3]. 




superconductivity, since the Bose-Einstein condensate has transition temperature that 
is inversely proportional to mass, and wavefunctions may not overlap. Such bipolarons 
are easily formed on lattices with triangular plaquettes in the presence of extremely 
large on-site Coulomb repulsion, and persist to large EPI. This conclusion is backed 
up by analytics in the anti-adiabatic approximation in the presence of large intersite 
Coulomb attraction. Another important conclusion is that the triplet-singlet exchange 
energy is of the first order in the hopping integral, and triplet bipolarons are heavier 
than singlets in certain lattice structures at variance with simple intuitive expectations. 
We summarise the types of lattices where light "crab" bipolarons may be formed in 
figure [121 contrasting with the traditional Holstein bipolarons (bottom) and describe 
the required effective interaction in figure [13] demonstrating the underlying physics of 
such bipolarons. 

Our CTQMC simulations lead us to believe that the following recipe is worth 
investigating to look for very high-temperature superconductivity: (a) The parent 
compound should be an ionic insulator with light ions to form high-frequency optical 
phonons, (b) The structure should be quasi two-dimensional to ensure poor screening 
of high-frequency c-axis polarized phonons, (c) A triangular lattice is required in 
combination with strong, on-site Coulomb repulsion to form the small superlight Crab 
bipolaron (d) Moderate carrier densities are required to keep the system of small 
bipolarons close to the dilute regime. Many of these conditions are already met in 
the cuprates. 
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